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ABSTRACT 

Eigenmode analysis is one of the most promising methods of analyzing large data sets in ongoing and 
near-future galaxy surveys. In such analyses, a fast evaluation of the correlation matrix in arbitrary 
cosmological models is crucial. The observational effects, including peculiar velocity distortions in red- 
shift space, light-cone effects, selection effects, and effects of the complex shape of the survey geometry, 
should be taken into account in the analysis. In the framework of the linear theory of gravitational 
instability, we provide the methodology to quickly compute the correlation matrix. Our methods are 
not restricted to shallow redshift surveys, arbitrarily deep samples can be dealt with as well. Therefore, 
our methods are useful in constraining the geometry of the universe and the dark energy component, 
as well as the power spectrum of galaxies, since ongoing and near-future galaxy surveys probe the uni- 
verse at intermediate to deep redshifts, z ~ 0.2-5. In addition to the detailed methods to compute the 
correlation matrix in 3-dimensional redshift surveys, methods to calculate the matrix in 2-dimensional 
projected samples are also provided. Prospects of applying our methods to likelihood estimation of the 
cosmological parameters are discussed. 

Subject headings: cosmology: theory — galaxies: distances and redshifts — galaxy clustering — 
large-scale structure of universe — methods: statistical 

1. INTRODUCTION 



o 
o 



in 

CO 



Oh 



O 

The large-scale structure (LSS) of the universe offers invaluable information on cosmology. The clustering pattern of 
galaxies reflects what the primordial universe looks like and how the universe evolves. In contrast to the observations of 
the cosmic microwave background (CMB), in which only information from photons emitted at a fixed redshift Zd cc — 1090 
are observed, the observations of the LSS provide information from photons emitted at various redshifts z 5. Therefore, 
O ' while the CMB primarily tells us the state of the early universe, the LSS shows us that of the present universe and 
+3 , its recent evolution. Detailed analysis of the LSS is extremely important in order to obtain a consistent picture of the 
• evolving universe. 

The spatial distribution of galaxies has particular importance in cosmology. Traditionally the power spectrum, or 
^ . the correlation function is computed from galaxy distributions and is compared with theoretical predictions. Strictly 
J^j ' speaking, the power spectrum and the correlation function are properly evaluated only when the sample is homogeneous. 
rS , In a real-world observation, homogeneous sampling of galaxies is not possible. The sampling density varies with radial 
■ distance, and possibly with survey directions. The survey geometry may have a complex shape. The observed galaxies 
are inevitably on a light-cone, therefore evolutionary effect comes in when a survey has a certain depth. 

Therefore, estimating the power spectrum from a real data needs a method to approximately erase those sources of 
inhomogeneity. Among others, Feldman, Kaiser, & Peacock (1994) offers a standard method to estimate the power 
spectrum. The estimated power spectrum by such standard method is a convolution of the true power with the window 
function of the survey. This becomes a problem when the survey has a complex shape. 

In a seminal paper, Vogeley & Szalay (1996) introduced a novel method, i.e., the method of Karhunen-Loeve (KL) 
eigenmodes in the context of galaxy redshift surveys. This method can handle complex shape of the survey geometry as 
well as complex selection functions. In the context of CMB analysis, Bond (1995) introduced an equivalent method which 
is called signal-to-noise eigenmodes method. The KL method is now recongnized as one of the most promising methods 
to tackle large data sets of the cosmological surveys (Tegmark, Taylor, & Heavens 1997; Tegmark et al. 1998; Taylor, 
Ballinger, Heavens, & Tadros 2001), and has been successfully applied to the Las Campanas Redshift Survey (LCRS) 
(Matsubara, Szalay k, Landy 2000), the Sloan Digital Sky Survey (SDSS) (Szalay et al. 2003; Pope et al. 2004), the ROSAT- 
ESO Flux-Limited X-Ray (REFLEX) Galaxy Cluster Survey (Schuecker, Guzzo, Collins, & Bohringer 2002). A variant, 
simplified method, which is called pseudo-KL eigenmodes method is applied to the IRAS Point Source Catalog Redshift 
Survey (PSCz) (Hamilton, Tegmark, & Padmanabhan 2000), the 2dF Galaxy Redshift Survey (Tegmark, Hamilton, & 
Xu 2002), and SDSS (Tegmark et al. 2003). 

In applying the KL technique to cosmology, it is important to have a good method to theoretically compute the power 
spectrum of the KL modes. This is not a trivial task, since the observed distribution of the galaxies are distorted by 
observational effects including peculiar velocity distortions, light-cone effects, varying selection functions, and complex 
shape of the survey geometry. Among these, the light-cone effect has never been taken into account in previous applications 
of the KL method to redshift surveys. The light-cone effects consist of evolutionary effects and geometric effects. The 
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clustering of galaxies evolves with redshift through evolution in mass clustering, and evolution in the biasing of galaxies. 
The comoving distance is not proportional to the redshift and the apparent clustering pattern in redshift space is not the 
same as in comoving space. For shallow redshift surveys in which the redshift range is, e.g., ^ z < 0.2, the light-cone 
effect has small contributions. However, QSO surveys probe much deeper regions, and ongoing and near-future redshift 
surveys, like the Deep Extragalactic Evolutionary Probe (DEEP) survey (e.g., Davis et al. 2003), the SDSS Luminous 
Red Galaxy (LRG) survey (Eisenstein et al. 2001), and the Kilo-Aperture Optical Spectrograph (KAOS) project 1 probe 
the universe of redshifts 0.2 ^ z ^ 5. 

Therefore, it is crucial in cosmological analyses of large surveys to have a good methodology in applying the KL 
eigenmodes method to ongoing and near-future redshift surveys which are not necessarily shallow. In this paper, methods 
to efficiently evaluate the theoretical values of the correlation matrix, playing a central role in the KL analysis, are 
developed and described in detail. This paper is organized as follows. In §2 the Bayesian analysis of the large data set in 
the context of galaxy surveys is reviewed, and necessity of a fast method to compute the correlation matrix is explained. 
In §3 the central method of this paper is presented. Necessary quantities in describing the method arc introduced in §3.1. 
A fast method to compute the linear correlation matrix in redshift space is given in §3.2 and §3.3. The nonlinear Finger-of- 
God effect is discussed in §3.4. Fast methods to compute the linear correlation matrix in projected 2-dimensional samples 
are provided in §3.5. In §4 additional techniques in applying the present methods to the KL eigenmodes are explained. 
Finally, the presented methods are summarized and discussed in §5. 

2. THE EIGENMODE ANALYSIS AND IMPORTANCE OF THE CORRELATION MATRIX 

In this section we briefly review the eigenmode analysis originally introduced by Vogeley & Szalay (1996) in the context 
of the redshift surveys. The importance of accurately modeling the correlation matrix in redshift space and of developing 
a numerically fast algorithm to calculate the model matrices, which are the main subject of this work, will be clarified in 
this section. 



2.1. Bayesian Analyses of the Huge Data Sets 

One of the central motivation of the eigenmode analysis comes from the difficulty in direct Bayesian analysis of huge 
data sets of observed galaxy distributions. It is desirable that the cosmological models are discriminated by the galaxy 
distribution itself, not relying on the inverse problem of determining the underlying density power spectrum from observed 
galaxy distributions. While a cosmological model predicts the linear power spectrum, the dependence of the galaxy 
clustering in redshift surveys on the theoretical power spectrum is not simple. 

Nonlinearity in the gravitational evolution distorts the power spectrum on small scales 10 h~ 1 Mpc. Peculiar velocities 
anisotropically distort the power spectrum on all scales. Spatial variation of the mass-to- light ratio or the galaxy biasing 
distorts the power spectrum. In deep redshift surveys, the spatial curvature and the evolution of the Hubble parameter 
distort the apparent pattern of the galaxy clustering. These various distortions in the observed power spectrum also 
depend on the cosmological model. Even if one can invent and adopt some ansatz on such a complex inversion problem, 
fair amount of information will be lost in this inversion process. 

However, the reconstruction of the underlying power spectrum is not necessary to discriminate the cosmological models. 
In the Bayesian approach, the posterior probability of the set of model parameters = {8 a }, given a set of observational 
data D = {di} and additional prior information I, is estimated by Bayes's theorem: 

P(G\DJ)=P(0\I) P W®> I ) \ (2.1) 

The first factor is the Bayesian prior, i.e., the prior knowledge on the model parameters, and the denominator is the 
normalization constant. Evaluating the numerator is the crucial part of the Bayesian approach. One needs to repeatedly 
evaluate this factor to maximize the posterior probability in multi-dimensional parameter space. Therefore, the numerator 
should be numerically calculated as fast as possible. This part of the calculation is usually not easy, especially when the 
the data D consists of a large set of observational figures. In galaxy redshift surveys, the data D consists of the distribution 
of galaxies in redshift space. These days the scales of redshift surveys become more and more immense. In the 2dF Galaxy 
Redshift Survey (Collcss et al. 2001), redshifts of about 220,000 galaxies are observed in the sky coverage of 2,000 square 
degrees. In the Sloan Digital Sky Survey (SDSS) (York et al. 2000), redshifts of about 1,000,000 galaxies are planned to 
be observed over 10,000 square degrees of the sky. 

When the prior information on the model, /, has nothing to do with obtaining the data themselves, the numerator of 
equation (2.1) is equal to P(D\0). In addition, when the statistical distribution of the data D is given by the multivariate 
Gaussian distribution, this factor has the form, 



P(D\0) = ^exp 
^(2tt) n dctR(0) 



-i(d T -(d T ) ) R~\0) (d-(d) ) 



(2.2) 



where d is the data vector d — (di, 0I2, • • ■ , d;v) and R is the N x N correlation matrix given by 



R(0) = ((d - (d) ) (d (d) e f) e = (dd T ) e (d) e (d T ) & . (2.3) 





http://www.noao.edu/kaos/ 
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In the above notations, (• • ■)© indicates the ensemble average over the statistical realizations of the universe (see, e.g., 
Peebles 1980) for an assumed set of model parameters 0. The correlation matrix therefore explicitly depends on the 
model parameters. 

In the analysis of redshift surveys, the observed set of data is primarily a set of the positions of galaxies in redshift space. 
The density field is derived from this primary data, which is one of the most fundamental information of our universe in 
redshift surveys. Therefore, it is natural that the data D is identified as the number density field of galaxies in redshift 
space. For this purpose, the redshift space is pixelized to obtain a discretized set of galaxy counts. Thus, the set of the 
number counts di in a cell i is the data D in redshift surveys. The pixel size determines the resolution in the analysis, so 
it is desirable to adopt as many cells in redshift space as possible. However, too many cells cause trouble in calculating 
the inverse of the correlation matrix in equation (2.2). The caluculation of the inverse of the huge matrix requires a lot 
of computation time, which is primarily an TV 3 process and can be an N 2 - 8 process at most (Press et al. 1992). This is 
the bottle neck in the Bayesian analysis since the inversion should be repeatedly performed for every set of parameters. 

The pixelized raw data vector d = (d\, cfo, • • • , g?at) t consists of clustering signals and noises. Most of the process in 
inverting the correlation matrix in the above procedure might be devoted to dealing with noises. Since the information on 
the cosmological model is only contained in the signals, it is beneficial to reduce the dimension of the correlation matrix, 
keeping the maximal information on the clustering signals and minimizing contributions by the noise. One simple way to 
achieve this aim is to consider a linear projection of the raw data vector d into vector B of a smaller dimension: 

B = Pd, (2.4) 

where the projection matrix P has a dimension M x N, and the dimension of the vector B is M, where M < N. The 
correlation matrix of the reduced data vector is given by 

C{0) = ((B (B) ) (B (B) e f) e = PR(0)P T . (2.5) 

The linear transformation from d to B ensures the preservation of the Gaussianity of the distribution. 

A clever choice of the projection matrix P enables us to avoid the inversion of the originally huge correlation matrices 
without losing most of the information contained in signals. In such a case, the inversion of smaller correlation matrices for 
the reduced data suffices for the Bayesian analysis, and the computation time is reduced by a factor of about (M/N) 2 - 8 . 
In Appendix A, the construction of the projection matrix by the KL eigenmodes with maximal signal-to-noise ratio is 
comprehensively reviewed. In short, the projection matrix by KL eigenmodes picks the first M modes that have largest 
S/N ratio. 

The probability of having a reduced data set, Adduced = {B n }, given the set of model parameters, 0, is therefore given 
by 

P(D ieduced \0) = V{2n)M \ ctC{e) exp [-± K - <* T >e) C-\G) (B - <B> )] . (2.6) 

The above equation assumes that the procedure of obtaining data is independent of the assumed model parameters, so that 
the projection matrix should be fixed throughout analyzing and maximizing the likelihood function. However, the optimal 
choice of the projection matrix usually depends on the underlying correlation function, and thus on model parameters, as 
explicitly presented in Appendix A in the case of the KL-mode projection. Therefore, the projection matrix is initially 
constructed by a fiducial set of model parameters, which is fixed throughout the maximization of the likelihood function. 

While the reduced correlation matrix C constructed by KL eigenmodes is diagonalized for a fiducial set of model 
parameters, it should be noted that the correlation matrix C is not diagonal for other set of model parameters. Therefore, 
the inversion of the correlation matrix C still needs to be performed. The likelihood function is not biased by the initial 
choice of the fiducial model, because the fiducial parameters only determine the projection matrix P which is fixed 
throughout. If the initial choice of the fiducial model is far from the correct model, the projection matrix is not the 
optimal one, and therefore we obtain rather broader profile of the likelihood function around the maximum, because of 
picking low-S/N modes. One can iterate the likelihood maximization by choosing a fiducial model from the preceding 
estimation to acquire the narrower profile of the likelihood function. This is the procedure we apply in practical situations. 

The model correlation matrix C{0) of a reduced data set should be repeatedly calculated for various models, fixing the 
projection matrix P. One still needs to construct the N xN matrix R(0) as many times as C, because C(0) = PR(0)P T . 
It is only when this process is omitted that the parameter dependence on R is linear or polynomial. However, most of 
the important parameters are not linearly or polynomially dependent on the correlation matrix. On that ground, a good 
modeling and numerically fast construction of the correlation matrix R is essential. 

2.2. Bayesian Analysis of the Galaxy Distribution in Redshift Surveys 

Now we turn our attention to the application of the general method explained above to the analysis of the data in 
galaxy redshift surveys. The primary data set in redshift surveys is the galaxy distribution in redshift space. In order 
to directly deal with the galaxy distribution, first we pixelize the redshift space, and then build the data vector d from 
galaxy counts di in cells of the volume Vj. The expected counts Ni = (di) are given by an integral of the selection function 
n(s), i.e., expected number density field in a survey: 



J d 3 sKi(s)n(s), (2.7) 
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where Ki{s) is the kernel function of a cell i, which is centered at Sj, and normalized by / d 3 sKi(s) = Vi. For example, 
the top-hat kernel of smoothing radius R is given by Ki(s) = 0(R — \s — Si\) and Vi = A-kR 3 /i, where 0(x) is the 
Heaviside step function and is the center of the cell i. We require the knowledge of the selection function in a survey. 
The correlation matrix is given by the moments of the count-in-cells (Peebles 1980), and is related to the convolution of 
the two-point correlation function in redshift space, £^(s, s'). Including the shot noise term, the relation is 

Rij = J d 3 sd 3 s'K l (s)K 1 (s r )n(s)n(s r )C {s) (s,s')+ J d 3 ■ 8 K i (a)K j (a)n{a) + E ij: (2.8) 

where the first term corresponds to the signal, second term corresponds to the shot noise, and the last term, Eij, is the 
correlation matrix for other sources of noise such as magnitude errors or uncertainty in the luminosity function, which 
is assumed to be independent from the clustering signal for simplicity. One can reasonably assume that the selection 
function n(s) within each cell is approximately constant and can be replaced by = iVj/Vj, in which case the equation 
(2.8) reduces to 

Rij = NiNj$ + ^NW 3 K t3 + E^ , (2.9) 

(no sum over i, j) where 

= Vv I d 3 sd 3 s'Ki(s)Kj(s')^ {a,a') , (2.10) 

and 

Kij = — j= J d 3 sKi( S )Kj(s). (2.11) 

The first term in equation (2.9) contains the model parameters, and corresponds to the signal part of the correlation 
matrix, Sij. The second term corresponds to the shot noise. When the cells do not overlap, the matrix is diagonal: 

— J d 3 sKi 2 (s), (non-overlapping), (2-12) 

(no sum over i, j), where Sij is the Kronecker's delta. For the non-overlapping top-hat kernel, this equation is simplified 
as = Sij , and the correlation matrix is simply given by 

Rij = NiNj^f + N l 5 lj + E^, (non-overlapping, top-hat) (2-13) 

(no sum over i, j). 

Unlike the correlation function in real space, the correlation function in redshift space is neither isotropic nor homoge- 
neous. The correlation function in real space is a function of separations between two-points, £(\x\ — x 2 \). However, the 
redshift space is distorted from the real space, and consequently the correlation function is no longer the function of the 
separation only. 

The comoving distance to a galaxy is not a primary observable quantity and the redshift is the only estimate of the 
real distance from the Earth to the cosmologically distant galaxies. Since the redshift is not identical to the comoving 
distance, the redshift space is anisotropically distorted from the real space. There are two predominant sources for the 
rcdshift-spacc distortion. 

The first source of the distortion is the peculiar velocities of galaxies. The peculiar velocities along the line of sight 
increase or decrease the observed redshift by the Doppler effect. On nonlinear scales, the peculiar velocities are dominated 
by the random motion in a cluster potential. On linear scales, the coherent motion along the density gradient is predom- 
inant. The linear velocity distortions are relatively simple when the distant-observer (or plane-parallel) approximation is 
applied. Kaiser (1987) found the linear power spectrum in redshift space, p( s \k) is anisotropically distorted along the 
line of sight, and is given by the formula, 

pV{k) = (l+0rf,) 2 P(k), (2.14) 

where P{k) is the power spectrum in real space, Hk is the cosine of the wavevector k to the line of sight. The counterpart 
of this formula for the linear two-point correlation function, £^(r), which is given by the Fourier transform of equation 
(2.14), was derived by Hamilton (1992). 

The second source of the distortion is what is called Alcock-Paczyhski (AP, hereafter) effect. Alcock & Paczyhski (1979) 
pointed out that the spherical objects in deep redshift region appears to be distorted in redshift space because of the 
nonlinear relation between the redshift and the comoving distance, and of the nonlinear relation between the apparent 
angular size and the angular diameter distance. Suppose that there is a spherical object of diameter I in units of comoving 
length. The redshift difference of this object is Az = H(z)l/c, and the apparent angular size is A6 = l/D\(z), where 
H(z) is the redshift-dependent Hubble parameter and Da(z) is the comoving angular diameter distance [proper angular 
diameter distance is given by d\{z) — Da(z)/(1 + z)}. The ratio 

_ Da(z) , . 

zA9 cz/H(z) ( ' ' 

is unity for sufficiently small z. In deep redshift regions, however, the ratio significantly deviates from unity and the 
deviation pattern strongly depends on the density parameter J?mo and on the cosmological constant J?ao, or parameters 
of the dark energy. 
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The redshift-spacc distortions of the two kinds, i.e., the velocity distortions and the AP distortions, are interwoven in the 
real universe except for the shallow redshift surveys without the AP effect. Ballinger, Peacock & Heavens (1996) derived 
the AP effect on the linear power spectrum, and Matsubara & Suto (1996) derived the same effect on the linear two-point 
correlation function, both using the distant observer approximation. Modern redshift surveys havee survey geometries 
with wide opening angles. In these cases, the distant-observer approximation is inappropriate and the wide-angle effect 
is not negligible in order to handle the information contained in the survey. Szalay, Matsubara & Landy (1998) derived 
a formula of velocity distortions in the linear two-point correlation function with wide-angle effect. Matsubara (2000) 
derived the most general formula of the linear two-point correlation function in which the velocity distortions, the AP 
effect, and the wide-angle effect are all included in a unified formula. All the previous formulas for the linear two-point 
correlation function are limiting cases of the last formula. Our construction of the correlation matrix R, which is explained 
in the following section, takes full advantage of the last expression. 

3. CONSTRUCTING CORRELATION MATRICES IN REDSHIFT SPACE 
3.1. The General Formula of the Correlation Function in Redshift Space 

To obtain the correlation matrix of equation (2.9), it is necessary to have the analytic expression of the correlation 
function ^ s \s,s'), which is given by Matsubara (2000). For completeness, we reproduce the formula in this section. 
Before that, we need to introduce some basic notations. 

The unperturbed metric is given by the Robertson- Walker (RW) metric, 

ds 2 = -dt 2 + a 2 (t) [dx 2 + S K 2 {x) (d6 2 + sin 2 9<j) 2 )] , (3.1) 
where we employ the unit system with c = 1, and adopt a notation, 

sinh (V— Kx) 



S K (x) 



(K < 0), 



x > ( K = °)> (3.2) 

sin 



(K > 0), 



<K 

The comoving angular diameter distance is related to the last function by 

D A (z) = S K (x(z)) . (3.3) 
where z is the "RW redshift" in which the peculiar velocity effect is not included, i.e., the redshift of the unperturbed 
universe. The radial coordinate x is the comoving distance from the observer, and the relation to the RW redshift is given 

by 

dz' 



where H{z) — a/a is the time-dependent Hubble parameter. We consider the general case that the dark energy has an 
equation of state p(z) = w(z)p(z), in which case, 

H{z) = H 0] J(l + z^Qmo - (1 + z) 2 r?K0 + (1 + ^) 3 cxp ^3jf ^^jn^ , (3.5) 

where f^Mo is the density parameter of matter(s), J7qo is the density parameter of the dark energy, and J?ko = ^mo+^qo — 1 
is the curvature parameter. 

The linear growth factor D(z) is a growing solution of the ordinary differential equation, 

D + 2H(z)D-^H 2 (2mo{1 + z) 3 D = 0, (3.6) 

where the dot represents the differentiation with respect to the proper time t, and the proper time t is related to the RW 
redshift z by 

f°° dz' 

t = J z (i + z>)H(z>y (3J) 

It is useful to transform the equation (3.6) to the following set of equations (Matsubara & Szalay 2003), 

(3.8) 

where a = (1 + z) 1 is the scale factor of the universe, and 



H 2 

n M {z) - tt^ttC 1 + ^) 3 %o (3.io) 



H 2 



Q Q( Z ) = tStt cx P f 3 f ) ^ Q o (3.11) 
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arc the time-dependent density parameters of matter and dark energy, respectively. Determining both the growth factor. 
D(z), and the logarithmic derivative of the growth factor, f(z), is neccesary in the formula of the redshift-space distortions. 
Accordingly, the Runge-Kutta integrations of the set of equations (3.8) and (3.9) meet the requirement. 

The linear formula for the two-point correlation function in rcdshift space is presented in Appendix B. It has the form, 

2 2n 

Z<- a Hzi,ZjAj) = b i b j D i D j Y,^<i n Hxi,XjAj)Zjr ) ( x ij)> ( 3 - 12 ) 

where Zi and Zj are the redshifts of the two-points and % is the angle between them. The linear bias parameter b(z) as 
a function of the redshift is introduced. The distances Xi = x(z{), Xj = x(zj) arc comoving distances to the two points 
according to equation (3.4). The comoving separation Xij between the two points is calculated by equation (B29). The 

functions c|™^ are given by geometrical quantities and defined by equations (B30) (B38). The functions £^ are given by 
integrations of the power spectrum, and defined by equation (B27). Although the expression (3.12) is somewhat long, it is 

straightforwardly calculated once the single-variable functions f(z), D(z), x(z), b(z), ^ n \x) are calculated and tabulated 
beforehand for a fixed cosmological model. 

While the correlation function linearly depends on the power spectrum through the function £,[ n \x), the dependences 
on other cosmological parameters like ^mo, ^qo, w(z) are nonlinear. The bias parameter b(z) only appears as 2nd order 
polynomials 

3.2. Smoothing Integrations of the Correlation Function 

Evaluating the correlation matrix of the equation (2.9) requires the six-dimensional integration of the correlation function 
in general. Direct integration is not a good idea. It is extremely important to obtain the elements of the correlation matrix 
without involving time-consuming steps such as multi-dimensional numerical integrations, because we need to construct 
many huge correlation matrices, with various cosmological models. We need to develop a method to calculate them with 
algebraic and interpolating methods without numerical integration for each element. There is actually such a method 
with special choices of the shape of the smoothing kernel as shown below. 

The signal part of the correlation matrix of equation (2.10) is essentially the correlation function of the smoothed 
density contrast, 

4 S) = (4 s) (^a^)4 s) (^^)) (3-i3) 

where 

5 { ^{x,6,<j>) = J S K 2 (x')dx' sm6'd6'dcl>'W R (x,6,(f>;x' ,6' ,(f>')5&(x' ,6' ,</>') (3.14) 

is the smoothed density contrast in rcdshift space, which is convolved by an normalized kernel function W R . The 
normalization of the kernel is given by 



/ 



S K 2 {x')dx'sm6'd6'd<t>'W R {x,6,(f>;x',6',(f>') = 1. (3.15) 



The kernel function W R (x, 9, 4>; x' , 6' , <f>') is non-zero only when the comoving separation between (x, 9, <f) and (x 1 ', 9' , <j>') is 
comparable to, or less than R. Practically, R is much less than the Hubble scale, and the linear operator in the equation 
(B24) is commutable with the smoothing integration of equation (3.14), since the functions D(z), b(z), f(z), a(z), and 
Sk{x) are practically constants over individual cells, and the Laplacian is an Hermitian operator. Therefore, the equation 
(3.14) reduces to 

5 [ s\x{z),e,cl>)=D{z){b{z) + f{z) ^(A + 3K)-^S R {x(z),0,<f>), (3-16) 

where 

5 R (x, 0,4) = J S K 2 {x')dx' sin & d9' dcf>'W R (x, 9, 0; x',6', 4>')5{x' , 9', <j>'), (3.17) 

is the smoothed density contrast in real space. A complex shape of the smoothing kernel complicates further analytic 
treatment, and therefore slows down the numerical evaluation of the correlation matrix based on the analytic formula. 
It is desirable to choose the smoothing kernel as simple as possible. In this repect, the simplest choice of the smoothing 
function W R is the spherically symmetric function in comoving space, in which case, the kernel is expanded by complete 
set of orthonormal modes X (c.f., eq.[B10]), 

/k 2 dk 
— W(kR)X (k,x rcl ), (3.18) 

where x ro i is the relative separation in comoving units between (x, 9, <f)) and (x' , 9', <f)'). The equation (3.18) is inverted as 



W(kR) = 4tt / S K {x)dxW R {x)X {k 1 x). (3.19) 
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Since R is much less than the curvature scale, R <C \K\ 1 / 2 , the above equation accurately approximates to 

W{kR) = 4tt / x 2 dxW R (x)j {kx) = [ d 3 xW R {x)e~ ik - x (3.20) 



Thus, the function W(kR) is the Fourier transform of the kernel, or the spherical window function of the kernel. The 
popular choices arc the top-hat window function W(kR) = 3ji(kR)/(kR) = 3(sinfci? — kRcoakR)/(kR) 3 , and the 
Gaussian window function W(kR) = exp(— k 2 R 2 /2). Useful window functions are given in Appendix C, including the 
m-weight Epancchnikov kernel, which might also be useful. 

There is a useful addition theorem for X\ (e.g., Ratra & Peebles 1995; Matsubara 2000), 

oo 

X {k,x Icl ) = ^(2i + l)Xi(M)*l(*^)fi(cos0 ro i), (3.21) 

1=0 

where Xi is given by equation (B6), Pi(x) are the Legendre polynomials, 8 IC \ is the relative angle between the directions 
(8, (j)) and (#', </>'). The sum rule (3.21) holds even for a non-flat case K ^ as well as a flat case K = 0. The Legendre 
polynomials are also expanded by the spherical harmonics Y™ ', 

fl(cos6i re i) = 2^i E Yr(e^)Yr*(e',n (3.22) 

m=—l 

Combining equations (3.18), (3.21), and (3.22), the spherical kernel function Wr is explicitly represented by the arguments 
(x, 8, (j); x' , 8', 4>') to be substituted into equation (3.17). The density contrast 5(x,8,<fi) in real space is also expanded by 
normal modes as in equation (B21). In this way, the integration by (x' , 8', <j>') in equation (3.17) is explicitly performed 
by orthogonal relations of equations (B3), (B9), resulting in 

/1.2 JL, _ 
-^- T W{kR)5i m {k)Xi{k, x)Y l m {8, 0). (3.23) 

This equation is a familiar one in the Fourier-Bessel expansion of smoothed fields in the flat universe, and is a natural 
extension to the non-flat universe, in general. The expansion (3.23) has just the same form as in the unsmoothed field of 

equation (B21), with the replacement of 8i m (k) — > W(kR)Si m (k). Consequently, the similar derivation of the unsmoothed 
correlation function in redshift space applies and the result of the smoothed correlation function is just given by the 
replacement P{k) — ► W 2 {kR)P(k) in the final expression of the unsmoothed expression in Appendix B [eqs. (B27) 
(B38)]. As a result, the correlation matrix (3.13) in redshift space with spherical cells of a fixed smoothing radius R in 
comoving space is given by 

2 2n 

$ = bibjDiDj E ^ i^M^ (*«)> (3-24) 



n=0 1=0 



where 



cW f ,_ (-1)" fk 2 dk ^Wk) v _ tu ^ r2l 



Q '(x)= \ ' / — ^ iy \ X l (k,x)W 2 (kR)P(k), (3.25) 

M v ) s K (x) J 2tt 2 k(k 2 -4K) n y ' ' y ' w ' v ' 

the factor Ni(k) is given by equation (B7) and the functions c\ n \xi,Xj,6ij) are given by equations (B30)-(B38). These 
equations, i.e., equations (3.24), (3.25), (B30)-(B38) are sufficient to construct a correlation matrix. The functions 

are given by elementary functions and the numerical evaluations are very fast. Although the functions involve 
integrations, this one-dimensional function is a common function to obtain the correlation matrix for a fixed cosmological 
model. Once the functions are evaluated and tabulated, one can use appropriate interpolations without integrations 

cell by cell. Therefore, one can make the evaluation of £,\ n \xij) a very fast procedure. These form our main technique of 
the fast construction of the correlation matrix. 

The simplicity of constructing the correlation matrix is largely dependent of restricting ourselves in linear theory. On 
nonlinear scales, the shape of the power spectrum is evolving, and the redshift distortions are not properly described by 
predictions of linear theory. To ensure the nonlinearity does not have significant effects, one should carefully choose the 
KL modes which are described by linear theory. This is achieved by calculating Fourier spectra of KL modes in indivisual 
survey (see, e.g., Pope et al. 2004). 

In practice, the comoving separations Xij are usually not comparable to the curvature scale of the universe. Observa- 
tionally, |J7ko| is less than 0.04 (Bennett et al. 2003; Spergel et al. 2003), which means the curvature scale is approximately 
five times larger than the Hubble scale, |-ftT| -1 / 2 £ 5c/Ho. The clustering scale we observe is usually much smaller than 
the Hubble scale. In this case, we can safely put K — > in the equation (3.25), since <§C |_ftT| -1 / 2 . That is, 

r (*) - J k ^^W 2 (kR)P(k), (3.26) 

when x <C |if | -1 / 2 . Even in this approximation, the effect of curvature on the separation between the observer and the 
sample, Xi = x(zi), are still included through the coefficients cj™\ The terms with a factor \K\SK 2 (xij) in cj n ^'s can also 
be neglected when the above approximation x^ <C \K\~ X I 2 is applied. 
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In the shallow rcdshift surveys, where z% <C 1 for all cells, one can totally neglect the curvature effect not only in Xij, but 
also in x^. Consequently, the relative distance x^ is simply given by x^ 2 = x 2 + Xj 2 — 2xtXj cos 6^. The distances are just 
given by the redshifts, Xi = czi/H , which does not depend on any cosmological parameter when distances are measured 
in units of h~ 1 Mpc. The quantity 7^ defined in equation (B40) is reduced to 7^- = (xi — Xj cos6ij)/xij and the quantity 

Oij in equation (B41) is simply given by 9^ = 6^. In the expression of the coefficients c|"' ) of equations (B30)-(B38), 
one can replace (3i — f3j = [3q = (3(0), and the terms proportional to li^S^ (x^) are negligible. The function a(z) in 
the equation (B25) or (B26) is predominantly determined by the variation of the selection function < P(z) or n(z) in the 
sample, so that the factor A^ in the equation (B39) is given by 




din $(z) 



dlnz 



x^ dlnn(z) 



Xi dlnz 



(3.27) 



The resulting formula of the correlation function is equivalent to the expression derived by Szalay, Matsubara & Landy 
(1998). The simplicity of the shallow surveys is that the correlation matrix depends polynomially on (3q = /(0)/6(0), and 
that there arc no other dependences on the cosmological parameters, when the power spectrum is given. 

In the limit Xy <C Xj, Xj, the distant-oberver approximation is fulfilled, which drastically reduces the complexity of the 
formula (Matsubara 2000). In this limit, one can reasonably assume that the quantity H(z)D(z)f(z)n(z) in the equation 
(B26) is approximately the same at z% and Zj so that = in the equation (B39). In the same limit, the approximation 

"fij = 7r — 7-,-i, and cos 9ij = 1 arc followed. As a result, the coefficients c| n ^ of equations (B30)-(B38) survive only for 

4 0) = l + ^(A + &) + ^&. (3-28) 



„(!) 



P 2 (cos 7 ), (3.29) 



cf = 35^/3^4(0037), (3.30) 

where Pi are the Legendre polynomials, and 7 = 7^ is the angle between the line of sight and the separation x%j . Although 
the distant-observer approximation is simpler, it is not appropriate in wide-angle rcdshift surveys. The general formula has 
more terms than the distant-observer approximation. However, this does not bring any significant problem to numerical 
evaluations of the correlation matrices. 

As an illustration, the resulting smoothed correlation functions calculated by our methods are plotted in Figure CI. 
We assume the cosmological parameters are those of the concordance model, Omo = 0.3, Oko = 0, w = — 1, h = 0.7, 
/baryon = 0.15, b = 1, erg = 1, where the Hubble constant in units of lOOkm/s/Mpc, h, the baryon fraction /baryon, and 
the normalization of the mass fluctuations as are needed only for modeling the CDM + baryon power spectrum P(k) (we 
used the fitting formula by Eisenstein & Hu 1998). The smoothing radius is 15/i _1 Mpc. In the Figure, the contour lines 
of the correlation function around centers at redshifts z = 0.04, 0.1, 0.2, 0.3, 0.4, 0.5, 1, 1.5, 2 are plotted. The plots are in 
"z-space" where the observer is sitting at the origin (0,0). More precisely, the contour lines indicates the values of the 

(s) 

correlation matrix £>/, where Zi takes the above discretized values. The remaining variables Zj and 9ij correspond to the 
coordinates (zj sin 9ij, zj cos 9ij) in the Figure. 

At low redshifts, the deviations from the distant-observer approximation are apparent. Thoughtless application of the 
distant-observer approximation results in a fatal error. Althought one can still choose parameterization in low-redshift 
correlations to obtain better results, our method is free from this ambiguity. 

There are prominent bumps in the contours which makes the contours look peculiar. These bumps are attributed to the 
baryon oscillation in the power spectrum. In fact, the real-space correlation function £(x) has the corresponding bump 
on the scale of around 100/i _1 Mpc, as shown in Figure C2. The thick curves correspond to the correlation function of 
the model which is equivalent to that used in Figure CI. There is a sharp peak on the scale of lOOft-^Mpc. This peak 
does not exist when /baryon = 0, as also shown in the Figure. 

When the smoothing is applied, the sharpness of the peak is reduced, and this "baryon peak" will be more like "baryon 
shoulder" in the correlation function. This baryon shoulder introduces a scale on the correlation function, which can be 
used as a "standard ruler" to measure the geometry of the universe (Matsubara & Szalay 2003; Blake & Glazebrook 2003; 
Seo & Eisenstein 2003). The location of the baryon peak in £(x) is sensitive to the parameter f^Moh which determine 
the horizon size at the time of equality. The amplitude of the peak is sensitive to the parameter /baryon which determine 
the strength of the baryon oscillation. Because the geometry of the universe determines the apparent scale of the baryon 
shoulder, the observation of £y constrains the geometry, almost irrespectively of the value of /baryon (Matsubara & Szalay 
2002). 

3.3. Spherical Approximations of the Kernel 

In the construction of the correlation matrix in the previous subsection, we have assumed the sphericity of the kernel in 
comoving space. In reality, the comoving distance-redshift relation and the curvature of the universe is needed in order 
to define a cell in redshift space, that has spherical comoving shape. This means that we have to assume the parameters 
^moj ^Qo> an d w(z) to exactly set the spherical shape. If we assume a wrong geometry, the assumed comoving space is 
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elongated or squashed to the line-of-sight direction. As a result, we end up with spherioidal cells in comoving space. One 
of the important application of calculating the correlation matrix in this paper is to determine cosmological parameters 
by Bayesian analysis, so that the uncertainty in setting spherical cells above is not desirable. However, the calculation 
of the correlation matrix for spheroidal cells is more complicated than spherical cells. The higher-order multipoles of the 
spheroidal cells cause an additional convolution in the expression (3.23). The simplicity of the correlation matrix is not 
preserved in that case, and thus the numerical computation of the correlation matrix will be significantly slowed down. 

However, such elongation or squashing does not significantly change the value of the correlation matrix. The elements of 
the correlation matrix for pairs of distant cells are virtually unchanged, since the correlation function is a smooth function 
and the elements of the correlation matrix for distant cells are virtually the same as a bare values of the correlation 
function without smoothing. The elements of the correlation matrix for nearby cells are affected by a shape of the cells, 
although the effect is not expected to be large as long as the elongation or squashing is not large. 

Keeping the simplicity of the expression of the correlation matrix, we introduce an approximate method to calculate 
the correlation matrix using spherical cells even when the assumed geometry is not exactly true. In this method, first 
we adopt a best guess of the geometry with assumed set of parameters, f?M0i ^qo> an d w(z). The KL-mode projection 
in the Bayesian analysis explained in the previous section also needs a fiducial cosmological model, so that it might be 
natural to take the same parameters as for that fiducial model, although it is not necessary. Next we set the cells in 
observed redshift space to be spherical in the comoving space of the assumed geometry, and construct a data vector by 
galaxy counts in these cells. The assumed model to construct the data vector will be called the "map model" below. The 
cell radius R is defined in this map model. 

In the Bayesian analysis, the data vector is fixed throughout the parameter estimation where cosmological parameters 
are varied. The parameters in the map model are generally not the same as the various cosmological parameters to be 
tested in the analysis. Therefore, the shape of the cells in the constructed data vector is not exactly spherical in each 
cosmological model to be tested, and is generically distorted to spheroids. In our approximation, the spheroidal shape 
in each cosmological model is approximated to the spherical shape with the same volume as the actual spheroids. The 
volume of the spheroids are given by 

[ D A (z) l 2 H^){z)A^ 

Vhcroids W - |^(map) (z) J 3 > 

where ff( ma P)(z) and D^ nap) (z) are the Hubble parameter and the comoving angular diameter distance of the map 
model, and H(z) and Da(z) are the counterparts of the tested model in the parameter estimation. The spheroidal cell is 
approximated to the spherical cell of the effective radius, 

1/3 



Res(z) = 



D A (z) 


2/3 




_D { ™ p) (z)_ 




I H(z) \ 



R. (3.32) 



Since this kind of the shape correction is necessary only for nearby cells, and the effective radius is almost the same 
for nearby cells, we do not need to consider the smoothing integration of the cross correlation with different smoothing 
radius, in practice. Interpolation of the correlation matrices with several effective smoothing radii is sufficient when the 
map model and the tested model is not too different. In defining the map model, it is natural to use the same cosmological 
parameters as in constructing the projection matrix P. Therefore, when the analysis is iterated with respect to optimizing 
the projection matrix, it is preferable that the cosmological parameters of the map model are also iterated, even though 
it is not strictly necessary. 

3.4. Insignificance of the Finger-of-God Effects on Linear Scales 

While the linear velocity field induces only the coherent motion of galaxies, random velocities are also present in the 
nonlinear regime, and the clustering shape in redshift space is elongated in the direction of line of sight. This is the famous 
finger-of-God effect which is not described by the linear formula. How does the finger-of-God affect on the correlation 
matrix? Since the finger-of-God effect decreases with distances between centers of the cells, the maximal effect is on 
the identical cells, that is, on the variance a 2 (R,z) = . Consequently, we can estimate the maximal influence of the 
finger-of-God effect on the correlation matrix by the same effect on the variance. 

The finger-of-God effect can be analytically modeled by incoherent peculiar velocities on small scales (Peacock 1992; 
Peacock & Dodds 1994). We can safely use the distant-observer approximation in this case. When the incoherent motion 
along the line of sight is empirically fitted by an exponential distribution function, f(v) = cxp(— y/2\v\/a p )/(\/2a p ) (Davis 
& Peebles 1983), the power spectrum with the finger-of-God effect is given by (Park et al. 1994) 

P FOG (k,Li) = ^ f , (3.33) 

(1 + £:Vct f 2 /2) 2 ' 

where /x is the direction cosine of the wavevector k with respect to the line of sight, and o-p = cr p /H(z) is the rms 
of the displacement in redshift space by incoherent velocities. The element of the correlation matrix, including linear 
rcdshift-distortion, is therefore modeled by 
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where Xij — Xi — Xj, [i — k\\/k and k\\ is the line-of-sight component of the wavevector k. Adopting the polar coordinates 
k = k(sin 9 cos (j>, sin 9 sin <j), cos 9) in fc-space, and setting = (x± , 0, x\\ ), the 0- integration can be analytically performed 
to give 

& = / ' k ^rW\kR)P{k) J^cce(kx^)J (kx ±y /T=f) (yTJ^M^) ■ (3 ' 35) 

The dimensionality of the above integration will not be analytically reduced further for arbitrary separation Xij. The 
integrand is a strongly oscillating function, so that the two-dimensional integration will be a slow process. 

However, the correction of the finger-of-God effect is not nessesarily needed as we will see below. The finger-of-God 
effect erases the clustering on scales smaller than cr F . Obscrvationally, the velocity dispersion is approximately given by 
the figure <r p = 340km/s (Davis & Peebles 1983) for z = 0. This dispersion corresponds to cr F = 3.4 h~ 1 Mpc, which is 
smaller than the nonlinearity scale, R ~ 10 /i _1 Mpc for z — 0. Therefore, the finger-of-God effect is expected to be smaller 
than the nonlinear effect on the clustering on average. Higher redshift regions have lower a p so that the finger-of-God 
effect is less significant. 

To find more quantitative estimates of what is the finger-of-God effect on the correlation matrix, we compare the 
variance with and without the finger-of-God effect. The non-diagonal elements of the correlation matrix are less affected. 
The variance <J 2 (R) is given by substituting x± = X|| = in the equation (3.35), and the /x-integration can be analytically 
performed to give (Cole, Fisher, & Weinberg 1995) 



a\R) 



Arctan {ka-p / y/2) 




fc 4 cr F 4 /4 



W 2 {kR)P(k). (3.36) 



Taking the limit <7 F — > 0, a familar amplification factor 1 + 3/3/2 + /3 2 /5 appears (Kaiser 1987). Focusing only on the 
finger-of-God effect, we set (3 — for the moment. As an illustration, the variance in real space and in redshift space 
(with only the finger-of-God effect) are shown in Table 1. We take a power-law spectrum P(k) oc fc -12 which corresponds 
to observations in nonlinear regime. The normalization a{R = 8/i _1 Mpc) = 1, and the velocity dispersion <r p = 340km/s 
are adopted, which match the observations (Davis & Peebles 1983). The power spectrum on large scales deviates from the 
power-law spectrum. In Table 2 the results for the cold dark matter (CDM) spectrum with a fitting formula of Eisenstein 
& Hu (1998) are shown. The damping by the finger-of-God does not dominate the linear infall effect on larger scales. 
The linear infall effect enhances the variance by a factor 1 + 3/3/2 + /3 2 /5. When Qmo = 0.3 and 6=1, the enhancement 
factor is 37%, and when 6 = 2, the factor is 17%. These figures are always larger than the finger-of-God effect when 
R> lOfc^Mpc. 

The dynamical nonlinear effect on the variance is actually larger than the finger-of-God effect. In fact, the nonlinear 
correction to the variance with smoothing is of order ctl 2 , where cl 2 is the variance of the linear theory. For example, the 
nonlinear loop correction to the variance without smoothing is given by (<5 2 ) = ol 2 + 1.82ctl 4 (Scoccimarro & Frieman 
1996). From Table 1, we see that the order of nonlinear corrections are always larger than the order of finger-of-God 
corrections by about a factor of five. Even when the radius of the cell R is still on nonlinear scales, the KL-mode projection 
can only leave the linear regime. In which case, the nonlinear correction is not needed in the evaluation of the correlation 
matrix, neither is the finger-of-God effect. 

3.5. Correlation Matrices in Projected Samples on the Sky 

When only the position on the sky is catalogued and the redshifts of galaxies are not observed, the data vector is a set 
of galaxy counts on the 2-dimensional sky, projected along the line of sight with some redshift distribution function n(z). 
When the redshifts are totally unobserved, the selection function spreads over a wide redshift range. With photometric 
redshift data, the selection function can be chosen so that the selection is approximately uniform in some narrower redshift 
range and rapidly drops outside that range. 



Table 1 

The finger-of-God effect on the variance for the power-law model with n = —1.2. 



R(h~ L Mpc) 


5 


10 


15 
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2.3 
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0.32 


0.19 
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3.69 x 10"* 




1.7 


0.59 


0.30 


0.18 


9.1 x 10~ 2 


3.66 x 10~ 2 


difference 


28% 


12% 


6.4% 


4.0% 


2.0% 


0.8% 
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One of the advantage of the projected sample is the absence of redshift-distortions. The correlation matrix is given 
by just integrating the angular correlation function w(9). Handling the correlation matrix is thus simpler than the 3D 
sample. When the 2D smoothing kernel is circularly symmetric, the elements of the correlation matrix are represented 
by just a function of the distance between centers of cells. Therefore, the advantage of the circular cells is obvious. We 
consider an efficient method to obtain 2D correlation matrices for the circularly symmetric cells below. 

In projected samples on the sky, the angular correlation function w(9) is given by 

/>OC 

w{6)= / dz 1 dz 2 n(z 1 )n(z 2 )D(z 1 )D(z 2 )b(z 1 )b(z 2 )£[x(z 1 ,z 2 ,9)}, (3.37) 
Jo 

where x(zi,z 2 ,9) is the comoving distance between two points with redshifts z\ and z 2 , and with an opening angle 9. 
This function is given by a geometric relation, 

S K 2 [x(z 1 ,z 2 ,9)} ^ S K 2 [x( Zl )} + S K 2 [x(z 2 )} ^ 2C K [x(z 1 )}C K [x(z 2 )}S K [x(z 1 )}S K [x(z 2 )}cos9 

-KS^ixiz^SK 2 ^)} (1 + cos 2 9) , (3.38) 

where x(z) is the comoving distance-redshift relation of equation (3.4). The selection function n(z) is the normalized 
redshift distribution dN/dz of galaxies with 

poo 

dz n{z) = 1. (3.39) 



The volume factor is included in n(z) so that the relation to the selection function per comoving volume 'P(z) is given by 
H(z)n(z) oc Sk 2 [x(z)]<P(z). The correlation function in real space £(x) is given by 

. f k 2 dk sinfcx . , . 

tto = Jl^Hh&) P W> (3 ' 40) 

which is derived from equations (B21), (3.21) and (3.22) [the equivalent formula of £(x) for an open universe is given in 
Wilson (1983)]. 

The direct 2-dimcnsional integration of equation (3.37) is numerically straightforward and relatively fast, once the 
3-dimcnsional correlation function £(x) is tabulated and interpolated. In this integration, we should use a dense sampling 
for small separation, x(z\, z 2 ,9), for a fast calculation. For example, the integration variable z 2 for a fixed z\ can be 
transformed to y± = \ z 2 — z\ \ + 2z\ sm{9/2). According to the sign of z 2 — z\ we use y + and y_ respectively. Logarithmically 
uniform sampling in ^-integration can achieve the dense sampling for small separations. In this paper, we follow this 
strategy in numerically evaluating the equation (3.37). The result of the numerical integration for a concordance model, 
^mo = 0-3, f?Ko — 0, w = — 1, h = 0.7, /baryon = 0.15, b = 1, ag, = 1 is shown on Figure C3 by a solid line. The selection 
function is simply taken as a step-wise function, n(z) — const, for 0.5 < z < 1 and n(z) — otherwise. In the angular 
correlation function, the baryon peak in £(x) is almost erased, although there is a hint of a bump around 100-200 arcmin. 

The signal part of the correlation matrix is given by the smoothing integration of w{9): 

= j sm6d9d<f> J sin 9' d9'd<j>'Wi(9, <t>)Wj{p' , <t>')w{9 Te i), (3.41) 

where Wi(9,(j>) is the smoothing kernel of a cell i with a unit normalization, and # ro i is the angle between directions 
(9, (/>) and (9' , </>'). When the smoothing kernel does not have simple shape, or varies from cell to cell, the straightforward 
evaluation of the correlation matrix is given by directly integrating equations (3.37) and (3.41). Firstly, the correlation 
function £(x) is calculated and tabulated. Secondly, using the interpolation of the tabulated £(x), the two-dimensional 
integration of equation (3.37) yields the angular correlation function w(9) which is again tabulated. Finally, the four- 
dimensional integration of equation (3.41) provides the correlation matrix. 

However, the above procedure is not a fast process because the smoothing integration involves a four-dimensional 
integration cell by cell. In the philosophy that the numerical computation of the correlation matrix should be maximally 
reduced, it is better to choose the shape of the cells as in the three dimensional analysis of the redshift surveys. The 
simplest choice is the circularly symmetric kernel Wg B (9) where 9 is the relative angle from a center of a kernel, 9 S is the 
angular smoothing scale. The normalization of this kernel is given by 

/>oo 

2tt / sm9d9W Bs {0) = 1 (3.42) 



Wi 



Table 2 

The finger-of-God effect on the variance for the CDM model with F = 0.2. 
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difference 


22% 
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In this case, the elements of the correlation matrix Wij depend only on the angular separation of the cells, 0y, i.e., Wy = 
u(9ij). Once the function u{9) is calculated and tabulated, the correlation matrix is obtained by simple interpolations, 
which are very fast procedures. 

The straightforward representation of the function C(9) is provided by the standard multipole expansion: 

^) = Y,^ r Pi(^s9)\W l \ 2 C l , (3.43) 
1=0 

where C; is the angular power spectrum and Wi is the multipole of the circular window function 

Wi=2tt [ sm9d9We s (9)Pi{cos9), (3.44) 
Jo 

The angular power spectrum is related to the spatial power spectrum P(k) by 

f°° k 2 dk 

C l=^J o ^mk)\ 2 P(k), (3.45) 

where ^i(k) is defined by 

&i(k)= dzn{z)D(z)b(z)X l [k,x(z)}, (3.46) 
Jo 

and x(z) is given by equation (3.4). Numerical evaluation of this expression is not so fast, because integrations over k and 
z and sum over I are involved for oscillating functions, Xi and Pi. In addition, numerical evaluations of Xi(k, x) might be 
tricky for a non-flat universe. Thus, it is better to seek a more efficient method to numerically evaluate the function u{9). 

If the angle 9 is sufficiently small and the flat-sky approximation can be applied, the angular correlation function is 
given by Limber's equation (Peebles 1980) 

f°° f°° kdk 

w{9)= dzH{z)n 2 {z)D 2 {z)b 2 {z) — J [k9D A (z)]P(k), (3.47) 
Jo Jo 2lT 

where Da{z) = Sk[x{z)\ is the comoving angular diameter distance, and J v (x) is the Bessel function. This approximation 
is compared to the exact integration of equation (3.37) in Figure C3. As one can see, Limber's equation is reasonably 
accurate for 9 ^ 100 arcmin. ~ 1.7 degrees within 0.5% for the concordance model. 

The two-dimensional Fourier transform (in the flat-sky approximation) of the angular correlation function is given by 



w(l) 



[ d^ wi 9) = 2n r 9dO M W M 0) = f ^g^ffflp ( -L-) , (3.48) 
J Jo Jo d cA (z) \d cA (z)J 

where the completeness relation of the Bessel function (e.g., Arfkcn & Weber 2001), 

xdxJ v (ax)J v {bx) = -5{a - b) {v > -1/2) (3.49) 



/ 

Jo 



is used. The two-dimensional Fourier window function is given by 

W(19 S ) = J d 2 9e- a - 9 WeM = ^ J 9d6J {W)We s {9), (3.50) 

which should not be confused with the three-dimensional one in the previous subsection although we use the same notation. 
It is straightforward to show that w(l) ~ C; and W(19 S ) ~ Wi in the limit of Z 3> 1 because ^(cos^) ~ Jo(l9) in the same 
limit. The top-hat window function corresponds to W(19 S ) — 2Ji(W s )/(W s ). The smoothed angular correlation function 
u{9) is the convolution of the angular correlation w{9) with the smoothing kernel Wg s (6) for both ends, and consequently, 
using the convolution theorem, we have 

r°° idi 

u{9) = / J Q (W)W 2 (W s )w(l) 
Jo 27T 

r°° Ufjh 

dzH(z)n 2 (z)D 2 (z)b 2 (z) —J [k9d cA (z)]W 2 [k9 s d cA (z)\ P(k). (3.51) 
Jo 27T 

The numerical performance of this two-dimensional integration is reasonably fast because the integrand of the variable z 
is a smooth function. The integration by the variable k is fast as long as the effective spectral index of P(k) is negative 
on the relevant scales which is indicated by the angle 9. 

Limber's equation is not appropriate for wide-angle correlations. However, the smoothing angle 9 S is much smaller than 
the opening angle of the survey. In this case, we can apply flat coordinates within the smoothing cells, and the function 
u{9) is straightforwardly represented by 

u{9) = JokWkI^! J 9 2 d9 2 d(f> 2 WeM)WeM) 

x w |V(0 + 2 cos cj) 2 - 0i cos (j)i) 2 + (6> 2 sin (f> 2 - 9 1 sin </>i) 2 ] . (3.52) 
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This four-dimensional integration is not needed for a separation angle larger than the size of the cells, 9 3> 9 Sl because 
the smoothing effect is small in that case, so that u{9) = w(9) in the lowest (Oth) order. Taylor expansion of the equation 
(3.52) with respect to 9\j9 and 9 2 /9 yields higher-order corrections to this simple approximation of the Oth order. The 
smoothing kernel of the form, 

We e {e) = s - 2 W S {6/6 S ), (3.53) 

has a scaling with respect to the smoothing angle, where W s is a single- variable function, independent on smoothing angle 
9 S , and normalized by 

poo 

(3.54) 
(3.55) 



Defining scaled moments of the kernel, 



2tt / tdt W s {t) = 1. 



m k = 2ir I tdt t k W s (t), 
Jo 



the Taylor expansion of the equation (3.52) up to forth order is straightforwardly calculated to give 



u{6) = w{6) + ^ 



/'((?) + 



w'(6) 



+ 



m 4 + 2rri2 2 
32 



'(<?) + 



2w"'{9) w"{6) w'{9) 



9 2 



+ 



e 3 



+ 



(3.56) 



For a top-hat kernel, W$ e (0) = 0(9 S — 9)/tt9 s 2 1 W s (t) = &(l — t)/n and = l/(fc/2 + 1) so that the numerical coefficients 
of second-order and forth-order terms are 1/4 and 5/192, respectively. 

On Figure C4, the smoothed angular correlation functions u{9) in the various approximations given above are compared 
when a top-hat kernel is adopted with a smoothing angle 9 S = 10 arcmin. In the range of angles plotted here, the Limber's 
equation is accurate enough and the equation (3.51) gives practically exact result. The first-order approximation, in 
which the terms up to 9 S 2 in equation (3.56) are kept, well reproduce the exact result for 9 £ 1.69 s . The second-order 
approximation does not significantly improve the first-order approximation. At 6 = 29 s , fractional errors are 4%, 0.9%, 
and 0.3% for truncations of the Taylor series up to zeroth, first, and second orders, respectively. Therefore, one does not 
need to perform the four-dimensional numerical integration even for neighboring cells, as long as cells do not overlap. We 
only need to numerically evaluate the bare angular correlation function w(9). Since the angular correlation function is 
practically a smooth function, the numerical evaluations of derivatives of w(9) are stable enough once w(9) is obtained. 

Thus, in practice, only the self-correlation wu = w(0) should be evaluated in some way other than the equation 
(3.56). This should be done only once, so that the direct four-dimensional numerical integration of equation (3.52) is still 
reasonable. Alternatively, since the Limber's equation is appropriate for the self-correlation, one can use the equation 
(3.51) with 9 = 0, 

f' 00 /"OO J 77 

«(0)= / dzH{z)n 2 {z)D 2 {z)b 2 {z) / — W 2 [k9 s d cA (z)} P(k), (3.57) 
Jo Jo 27r 

which is two-dimensional integration and is faster to compute than the direct four-dimensional integration. 

Obtaining the one-dimensional function u(9) by using appropriate methods explained above, the correlation matrix Rij 
including the noise term is immediately constructed by (c.f., eq.[2.9]) 

Rij = NiNjutfij) + ^NJ^K,, + En, (3.58) 

where 9ij is the separation angle between cells i and j, Eij is the correlation matrix for other sources of noise except the 
shot noise. The quantities Ni, Kij are defined by 



Ni= I sm.6ded4>Ki{e,4))n{e,(t)), 
1 



Ki 



Si Si 



sin 6d6d4>Ki{6, <j>) Kj (9, <j>), 



(3.59) 
(3.60) 



where n(9, <p) is the mean number density of the 2D sample, which can depend on the direction because of possibly 
inhomogeous sampling, Ki(6, (f>) is the smoothing kernel of the cell i with an arbitrary normalization, 

Si = J sin 9d9d<j>Ki(6, 0). (3.61) 

The normalized kernel Wi{0, (f) in equation (3.41) is given by Wi(9, <p) = Ki(9, <fi)/Si (no sum over i). 

4. APPLICATION TO THE LIKELIHOOD ANALYSIS 

4.1. Cosmological Parameters of Simple Dependendence on the Correlation Matrices 

In the likelihood analysis, one needs to repeatedly calculate the correlation matrices for models with various parameters. 
The methods to calculate the correlation matrices described above are already fast. However, when the dependence on 
some cosmological parameter is linear or polynomial, the repeated compression of the correlation matrix R into a reduced 
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matrix C of equation (2.5) is not necessarily needed in changing that parameter. Generally, if the dependence of a 
cosmological parameter 9 on the correlation matrix has the form, 

R = RW + R^9 + R^9 2 + ■■■, (4.1) 

where R^ 's do not depend on 9, then the reduced correlation matrix of equation (2.5) is given by 

C = C* (0) + C {1) 9 + C {2) 9 2 + ■■■, (4.2) 

where 

Cr(fc) = prWpt, (4.3) 

and P is the projection matrix of equation (2.4). Therefore, when the partial matrices are calculated fixing other 
parameters, one can omit the projection of the huge matrix in changing the particular parameter 9. Instead, the reduced 
correlation matrix C is obtained simply by a summation of equation (4.3). 

The obvious example of the cosmological parameter which has the polynomial dependence on correlation matrices is 
the normalization parameter of the mass power spectrum, cr 8 2 , or A s . The signal part of the correlation matrix, £y, or 
Wij is always proportional to this parameter. The correlation matrix of a rcdshift survey is written as 

Rii = ^NiN&j + ^/NjfjKij + Eij (4.4) 

where £y is the signal part of the correlation matrix of equation (3.24) with the normalization <r$ = 1. 

The dependence on the parameter f3 at a fixed redshift is polynomial. To see this, the expression of the correlation 
matrix £y of equations (3.24) and (B30)-(B38) is re-arranged according to the dependences on and (if 

& = bibjDiDj + + + (4.5) 
= DiDj fabd? + fibj^ + fM$ + fifjg?) (4-6) 

where ^ (k — 0,1,2) do not depend on or on b{, bj. The matrix ^ff is not symmetric in general. It is 

straightforward to obtain ^\ For convenience, below we give the explicit representation in the case of X{j <C l-ft'l -1 / 2 , 
which is fulfilled in practice, since the curvature scale is observationally more than five times larger than the Hubble scale 
as discussed just above equation (3.26). In this case we can put \K\Sk — 0, Cx(xij) = 1 and Oij = 9^ in the 

equations (B30)-(B38). Therefore, 

4 X) = + A v coamZPixij) + (cos 2 lt3 - Vj &\ Xij ), (4.8) 

4 2) = ^(l + 2cos 2 6 ij )$\x ij )-±A ij A ji cose ij £\x ij ) 

+ i [Aij (cos 7^ - 2 cos jji cos % ) + A,-* (cos 7^ - 2 cos 7^ cos % )] ^ (ly ) 



2 4 

- + - cos 2 9ij — (cos 2 7y + cos 2 jji) + 4 cos 7^ cos 7^ cos 9i 



£.2 \ x ij) 



+ A lJ A li (cos 7^ cos 7^ + icos%^j ^(xij) 

H — [Ay (5 cos "fij cos 2 7jj — cos "fij + 2 cos 7^ cos 0jj ) 
5 



1 

+ 7 



+ (5 cos 7jj cos 2 7^ — cos 7^ + 2 cos 7^- cos fly ) ] £3 (xy ) 
1 2 

- + - cos 2 9^ - (cos 2 7 i; ,- + cos 2 Tjj) + 4 cos 7^ cos 7,1 cos % 



+ 7 cos 2 7ij cos 2 7ji ] Ci 2) (sy )• (4.9) 

The quantity £^ corresponds to the isotropic component of the correlations, since it depends only on Xij. The quantities 

^ff and (ff are relevant to the distortions by the peculiar velocity field. 

When a fixed redshift z p i vo t is arbitrarily chosen, the correlation matrix is explicitly polynomial with respect to the 
parameter /3 piv = f3(z pivot ): 

in = {bpi v D piv ) 2 bibjDiDj + /3 piv faff + + (3 piv %P^ff} , (4.10) 

where b plv = 6(z p i vo t) and D plv = D(z plvot ) are the bias and the growth factor at a pivot redshift. The relative evolutions 
with respect to z p i V ot are represented by quantities b t = b(zi)/b(z pivot ), D { = D(zi)/ D(z pivot ), and % = (3(zi)/ f3(z pivot ). 
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When the redshift range of the survey is not significantly large, the relative evolutions are small, and are not strongly 
varying functions of redshift. 

It is not hard to imagine that the typical bias 6 p i vo t and the overall normalization of the mass power spectrum at z p i vot , 
o"8. P iv = -DpivCs, contribute to the correlation matrix in similar way, and these parameters degenerate in estimating the 
likelihood function, although the distortion terms break this degeneracy to some extent. Since it is not advantageous to 
deal with almost degenerated parameters separately, it might be preferable to consider the normalization of the galaxy 
power spectrum at z p i vot , 0"8g,piv = & p i v 08,piv = b p i v D p i v as, as an independent parameter instead of the unobservable 
a$. We introduce the normalized correlations, £^ (k = 0,1,2), which are calculated from the power spectrum with a 
normalization er 8 = 1, so that £^ = a 8 2 s}j\ The correlation matrix of equation (4.10) is then reduced to 

Rij = crsg^NiNjbibj DiDj + /3 P iv (M? + MS) + fav%M?] + + E ih (4.11) 

(no sum over i, j). Thus the normalization of the galaxy power spectrum at the pivot redshift, <78 g . p iv, and the redshift- 
distortion parameter at the pivot redshift, /3 p ; v , are independent parameters both of which polynomially contribute to the 
correlation matrix. The parameters a 8 and 6 p i v are dependent parameters through the relations 

/ (Zpivot ) D (z p i vot ) 

& P iv = %^- (4-13) 

Ppiv 

In similar ways, one can arbitrary choose two independent parameters out of cr 8 , 6 piv , /3 piv , and <78g,pivi depending on 
individual analysis. 

Each term in the equation (4.11) is separately projected to the reduced matrices as explained above. As detailed in 
Appendix A, the projection matrix P whitens the noise correlation matrix, and the reduced correlation matrix C = PRP T 
is given by 

C' nm — C8g.pi v [Cnm /^pivCnm Ppiv C^m^ ^nrrn (4-14) 

where 

c { °l = J2 r,j\,, J x,x J h,i b i>.i>A:: , (4-15) 

C« = E PniP^NiNjDiDj (f,hX,: + JMf ) - (4-16) 

- E PmPmjNtNjDiDjJ&f . (4.17) 

Once the projections of the matrices in the equations (4.15)-(4.17) are obtained for a given set of parameters except 
C8g,piv and /3piv, the reduced correlation matrix C nm is immediately derived by the equation (4.14) for any values of 
csg.piv and /3piv without projecting the huge matrix again. Other cosmological parameters, including J?mo, ^qo, and 

parameters which parametrize the evolution of the bias, b(z), and of the equation of state for the dark energy, w(z), are 
non-polynomially dependent on the correlation matrix Rij. The projection operation should be performed for each set of 
these non-polynomial parameters. 

In shallow redshift surveys, however, the evolutionary effects of the cosmological quantities above are small. In this 
case, the natural choice of the pivot redshift is z p i vo t = so that <78 g . p i v = &8g and /3 p i v = Po- One can also approximately 
set D = 1, f = 1, b = 1, and H(z) = Hq in the formula of the correlation matrix (Ho does not explicitly affect 
the correlation matrix in redshift space for a given power spectrum, only through the power spectrum). The only non- 
polynomial parameters are those which parametrize the power spectrum, since, with the good approximation, the distances 
are proportional to the redshifts and the geometry is flat. Therefore, in this case, the elements Qj are fully determined 
where any dependence on cosmological parameters is only through the power spectrum. The representation of equations 
(4.15)-(4.17) are simplified as 

C^ = Y,PniPm j N i N j ^\ (4.18) 
C£l = £ PniPmjNiNj (g) + g}>) , (4.19) 

hj 

C^l = Y J PniP mj N i N j ^\ (4.20) 

hi 

which are only dependent on the shape of the power spectrum. 
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4.2. Band Power Estimation of the Power Spectrum 

The correlation matrix is always linearly dependent on the power spectrum as long as the linear regime is considered. 
Therefore, decomposing the power spectrum into band powers can be used to straightforwardly estimate the power 
spectrum, without parameterizing the shape of the power spectrum (Tegmark, Hamilton, & Xu 2002). Below we briefly 
explain the band power estimation within the context of our method. 

The power spectrum is decomposed into band powers as 

rip 

P(k) = J2A p P p (k), (4.21) 

P =i 

where P p {k) is the band-power base function which has support near the wavenumber k p , and A p is the power of the band 
p. The choice of the band-power base function is not unique, as long as the power spectrum is parameterized by linear 
parameters A p . The simplest choice is the piecewise constant functions, 

Pp(k) = 0(k p+1 -k)0(k-k p ). (4.22) 
Another choice is the piecewise linear functions, 

p p(k ) = k - Y 1 0(k p - k)0(k - fc p _i) + t P+1 ~ * 0(k p+1 - k)0(k - kp). (4.23) 

K p Kp—l K p+l K p 

One can construct more complex bases by optimizing the band powers (e.g., see Hamilton 1997; Hamilton & Tegmark 
2000). Each band contributes linearly to the reduced correlation matrix C nm . We denote the contribution of each band 

p to the matrices ^ {a = 0, 1, 2) of equations (4.7)-(4.9) by ti^' 1 ^ such that 

Hp 

P =i 

Then, analogously to the equations (4.14)^(4.17), the reduced correlation matrix is given by 

2 ip 

C nm = E Wpi^Ct? + <*nm , (4.25) 

a=0p=l 

where 

Cf m P) = E PniPrnjNiNjWADj^ , (4.26) 

C n ]f> = E PmPmjNiNjDiDj (ihX,^ + fM^) , (4-27) 

Cil P) - E PniPmjN.NjDiDjfJj^ ■ ( 4 - 28 ) 

■»(a,p) 



The matrices Cnm still depend on the cosmological parameters because of the evolutionary effects. However, in case 
of shallow surveys, they become independent on any parameters at all. Once the matrices Cnm^ are obtained, the 
dependencies of the correlation matrix on parameters a p and (3q are explicitly polynomial. 



4.3. Quadratic Estimator 

When the number of parameters to be estimated is large, it is not realistic to calculate the likelihood function on every 
mesh point in the parameter space. When the behavior of the likelihood function is completely unknown, there are no safe 
methods to maximize the likelihood in parameter space. In the present application in cosmology, the correlation matrix 
smoothly varies with cosmological parameters in general. In this case, the likelihood function near a peak is approximated 
as a Gaussian. Therefore, as long as we approximately know the true cosmological parameters, one can reach the location 
of maximun likelihood by a few iterative steps using the Ncwton-Raphson method. This approach is taken by Bond, Jaffe, 
& Knox (1998) (hereafter, BJK) in estimating the power spectrum of the cosmic microwave background (CMB). The 
same method is also applied to the analysis of the 2-dimensional galaxy survey (Hutcrcr, Knox & Nichol 2001; Efstathiou 
& Moody 2001). The purpose of this section is to give an aspect of applying BJK's method to the present analysis. 

Following BJK, the likelihood function of equation (2.6), £({9 a }) = P(-D r cduccd|<9) is approximated by a Gaussian: 

In C ({9 a + S6 a }) = In C ({6 a }) + £ ^66 a + WW 50 * 50 ? ^ 

a a a, 13 a 

In this approximation, the parameters which maximize the likelihood is directly solved as 

d0 a d0p) ~W 



W » = -E bos: w (430) 
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where first- and second-derivatives of the log likelihood function is explicitly calculated from equation (2.6), resulting in 



d6 a 2 
d 2 ln£ _ 1 

d8 a de fj ~ 2 r 



-c- 1 c a c- 1 c. 



(4.31) 



(4.32) 



where C is the correlation matrix and standard notations, C Q = dC/d9 a , etc., are employed. Instead of intensive 
calculation of the curvature matrix of equation (4.32), the second-order derivatives are replaced by an expectation value, 
which is equivalent to the Fisher matrix: 

d 2 ln£ \ 



1. 



v d6 a d6p 

In this way, we obtain the BJK's quadratic estimator, 

1 



= 2 E ^"'U Tr [( BBT - C ) (C-'CpC^) 



(4.33) 



(4.34) 



which proves to be very useful. Only few iterations are needed to reach the local maxima of the likelihood function when 
the band power of the CMB spectrum is estimated by this method (BJK). Since the projection matrix P of the projected 
correlation matrix C = PRP T of equation (2.5) is fixed throughout the likelihood maximization, the derivatives C. a are 
straightforwardly calculated by numerically differentiating the correlation matrix R: C a = Pi? jQ P T . For the cosmological 
parameters with polynomial influence on the correlation matrix, as considered in the section 4.1, the derivatives C Q are 
obviously given without any numerical differentiation. 



5. CONCLUSIONS AND DISCUSSION 

Technical details of the methods in constructing the correlation matrix in redshift space, which are essential in direct 
likelihood analysis of the large-scale structure, have been presented. Making use of the most general analytic formula 
of the linear two-point correlation function in redshift space given by Matsubara (2000), a fast procedure to construct 
correlation matrices in redshift space has been explicitly provided. This procedure can be used for deep redshift surveys, 
as well as shallow surveys. The finger-of-God effects are not significant as long as the nonlinear modes are excluded by 
truncating the KL-modes. We have also given fast methods to produce correlation matrices in projected samples on the 
sky in which spatial and spherical curvature effects are taken into account. The Limber's equation is not appropriate 
for large separations. Some parameters polynomially depend on the correlation matrix. Those parameters are evidently 
simple to analyze. When the number of parameters are large, the quadratic estimator is a promising way to maximize 
the likelihood function. 

The purpose of this paper is to describe technical details of computationally quick construction of correlation matrices. 
Parts of the present methods have been successfully applied to the real data, including the Las Campanas Redshift Survey 
and SDSS (Matsubara, Szalay & Landy 2000; Szalay et al. 2003; Pope et al. 2004). So far the application has been 
restricted to the shallow samples and the projected samples of galaxies. One of the final targets of the present methods 
is the application to the large-scale redshift surveys, which are deep and wide enough in redshift space. In this respect, 
our methods are effective in the analysis of the ongoing SDSS luminous red galaxy (LRG) sample (Eisenstein et al. 2001), 
in which ~ 10 5 galaxies are catalogued with redshift range of z ~ 0.2-0.5 over ~ 10 4 square degrees on the sky. Since 
this sample is dominated by the shot-noise on smaller scales, the KL transform is essential to maximally extract the 
cosmological information. The rcdshifts in the LRG sample are not very shallow as in the main galaxy sample, in which 
z ?S 0.2, the evolutionary effects on the clustering is detectable, and therefore geometry of the universe, nature of dark 
energy, etc., can be finely constrained because of the cosmological redshift-space distortion effects (Matsubara & Szalay 
2001, 2002, 2003; Blake & Glazebrook 2003; Seo & Eisenstein 2003). In next-generation redshift surveys beyond the 
SDSS, the present methods will provide a unique technique in cosmological analyses of the survey data. In large and deep 
redshift surveys of the future, the linear regime will be more focused on, and the shot noise will be more severe than the 
present-day surveys. The advantage of the methods developed in this paper is much greater in those future surveys than 
the past shallow surveys. 
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APPENDIX 



CONSTRUCTION OF THE PROJECTION MATRIX BY THE KL EIGENMODES 

In this appendix, the construction of the projection matrix P of equation (2.4) which reduces the dimensionality of the 
data space is reviewed, based on the method of Vogeley & Szalay (1996). 
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To capture maximum signals with minimum noise, one needs to keep the modes with the high S/N ratios and discard 
the modes with the low S/N ratios. For this purpose, we decompose the correlation matrix R into the signal part S plus 
the noise part N, R = S + N, assuming signals and noises are mutually uncorrelated. To normalize the noise correlations, 
we consider a linear transformation, 

d! = Qd, (Al) 

where Q is a non-degenerate N x N matrix, det Q ^ 0, which is not neccesarily a symmetric matrix. The transformation 
Q is chosen so that the noise correlation matrix in the new data vector d! is transformed to an identity matrix: 

N' = QNQ T = I, (A2) 
This transformation is referred to by prewhitening (Vogeley & Szalay 1996). 

However, the choice of the prewhitening matrix Q is not unique when the noise correlation matrix is not diagonal. For 
example, one can first diagonalize the noise correlation matrix N and then multiply a weighted diagonal matrix: 

Q = Q*Qi, (A3) 
where Q\ is an orthogonal matrix which diagonalize the noise matrix N: 

Q\NQi T = diag.(cri 2 , cr 2 2 , . . . <jn 2 ), (A4) 
and Q2 is a diagonal matrix weighted by the inverse of the square root of the noise eigenvalues: 

Q 2 =dmg.{a 1 -\a 2 -\... ( j N - 1 ). (A5) 
This choice is not the only possibility. Another example of the transformation Q is given by the Cholesky decomposition 
of the noise matrix, N = LL T , where L is a lower triangular matrix. The choice Q = L^ 1 also satisfies the equation (A2). 
Unless the noise matrix has diagonal form in first place, the matrix Q2Q1 is not a triangular matrix and the matrix 
is a lower triangular matrix so that Q2Q1 7^ £ _1 in general. 

Once the noise matrix is normalized to an identity matrix by a certain choice of Q, the correlation matrix of equation 
(2.3) is transformed to 

R' = QSQ T + I. (A6) 

Now, to obtain the statistically orthogonal set of data vectors, the prewhitened data is rotated so that the correlation 
matrix R' is diagonalized, solving the eigenvalue equation, 

R'* n = \',M'„. ^ (A7) 

Since the diagonalization of the matrix R' keeps the noise term as the identity matrix, the noise contribution to an 
eigenvalue A' is always unity. Therefore, the larger the eigenvalue is, the larger the signal-to-noise ratio of the corresponding 
eigenmode is. Sorting the eigenvalues in decreasing order, 

A' 1 >A 2 >--->A' M >--->A' JV , (A8) 
and discarding the modes with low eigenvalues of A' M+1 , . . . , A^y, the natural choice of the projected data is 

B n = *ld! = *lQd (A9) 
where n < M, and the eigenvectors satisfy the orthonormality, 

= Snm . (A10) 

Employing the projection matrix of the whitened data, 

P' = (# u # 2 ,...,# M f, (All) 

the optimal projection matrix P is given by 

P = P'Q (A12) 

The eigenvalue X' n corresponds to "the signal-to-noise ratio plus unity" for each mode, since one can notice that the 
matrix Q is the "square root" of the noise matrix (c.f., eq.[A2]). Therefore, the quantities A„ = \' n — 1 correspond to the 
"signal-to-noise" eigenvalues. This is more clearly shown by the fact that the eigenvalue equation (A3) is equivalent to 
the generalized eigenvalue equation 

S$ n = \ n N& n , (A13) 

where ^„ = Q T ^ n is the signal-to-noise eigenmode, which is introduced by Bond (1995) in the analysis of the CMB data. 
The projection matrix P is simply given by 

P= (3>i, # 2 ,..., #m) T - (A14) 
In this representation, it is explicitely seen that the choice of the projection matrix P is independent on the choice of the 
prewhitening matrix Q. The signal-to-noise eigenvectors # are no longer orthogonal, but are normalized by 

&lN& m = S nm . (A15) 
A schematic explanation of obtaining the S/N eigenmodes is given in Figure C5. The solid ellipses represent the signal 
correlations and the dashed ellipses represent the noise correlations in two-dimensional data space. The S/N eigenmodes, 
or the KL modes, the directions of which is indicated by black arrows, determines the direction of the projection to reduce 
the dimension of the data space, retaining the maximal signals-to-noise ratio. 

In the above procedure, the correlation matrix of the reduced data B is a diagonal matrix since it is obtained by 
diagonalization of the correlation matrix of the prewhitened data d! . Therefore, the reduced data is statistically orthogonal: 

(A16) 

where \' n = ^R'^n = #^i?#„ = A„ + 1 = ((B n ) 2 ) - (B n ) 2 . The above orthogonality is achieved only when the true 
cosmological model of the correlation matrix R is known. In practice, the true parameters of the cosmological model 
are the object of study. Therefore the orthogonality of C nm is only approximate and is never assumed in our likelihood 
analysis. 
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THE LINEAR FORMULA OF THE CORRELATION FUNCTION IN REDSHIFT SPACE WITH THE CURVATURE EFFECTS 



In this appendix, we summarize the most general expression for the linear correlation function derived by Matsubara 
(2000). The high-redshift effects, the wide-angle effects, the peculiear velocity effects, and the selection effects are all 
included. 

To describe the formula, we need the orthonormal modes of the Laplacian in the spatial section of the RW metric. That 
is, the orthonormal modes of the equation 

(A + k 2 - K)Z = 0, (Bl) 
prove to be useful. In the spatial section of the RW metric of equation (3.1), the above equation reduces to 



S K \x) 



d 



dx 



dZ 



1tASk\x)—\ + - 



dx 



1 d 



sin 6 d6 



n dz , 



1 d 2 Z 



sin 



+ (k 2 - K)Z = 0. 



(B2) 



Separating the variable Z into a radial part and an angular part, the orthonormal modes have the form Xi(k, x)Y l m (9, <p). 
where Y™ are the spherical harmonics. The orthonormality and the completeness of the spherical harmonics are useful: 



/ 



sin 9d0d^Y t m * (9, <\>)Y^ {0,<f>) = 6 W % 



1=0 m=-l 

The radial part Xi satisfies the differential equation 



5(9 - 9')5(cj ) - </>') 
sa\9 



1 



d 



S K 2 (x) dx 



S K \x) 



dx 



1(1 + 1) 



S K 2 (x) 



x, = o. 



(B3) 
(B4) 

(B5) 



1 /2 

Putting / = Sk Xi and dz — Sxdx, this equation reduces to the associated Legendre differential equation. The 
solutions of this equation are given by the conical function, the Bessel function, and the toroidal function, for negative, 
zero, and positive curvatures, respectively (Harrison 1967; Wilson 1983; Matsubara 2000). These solutions which are 
regular at the origin x = are explicitly represented by 

S K \x) f 1 — ' 



X t (k,x) = 



d 



y/NJk) \S K (x)dxJ \S K (x) 



sinkx 



where 



Ni(k) = Y[(k 2 .fK) = k 2 (k 2 - K)(k 2 - AK) ■ ■ ■ (k 2 - l 2 K), 

3=0 



(B6) 



(B7) 



arc the normalization constants. For x — > 0, Xi behaves like x l with appropriate constants, and thus Xi(k,0) = <5;o 
(Harrison 1967). For the flat universe, K = 0, the equation (B6) is simply given by the spherical Bessel function, 
Xi(k,x) = (—l) l ji(kx). For K < 0, k takes all positive values. In the case of positive curvature, K > 0, the 3-space is 
periodic and only discrete values k = (I + 1)VK, (I + 2)^/K, (I + 3)\/K, . . . are allowed. In the latter case, the functions 
(B6) are reduced to the Gegenbauer polynomials. The following recursion relations hold (Matsubara 2000): 



y/tf-PKXt-ifax) + (21 + \)^-\xi{k,x) + ^k 2 -(l + l) 2 KX l+1 (k,x) = 



S K (x) 



The orthonormal and completeness relations are 



f 

Jo 



4tt / dxS K 2 (x)Xi(k,x)Xi(k\x) 
Jo 



2n 2 S(k-k') 



k 2 



k 2 dk 
~2^ 



X l (k,x)X l (k,x') 



S(x-x') 
AnS K 2 (x) 



for K < 0, and 



/>oo 

4?r / dx S K 2 (x)Xi(k n ,x)Xi(k n ,,x) = 
Jo 



2tt 2 cL 



2tt 2 



"Y] k n 2 Xi(k n ,x)Xi(k n ,x') 



k 2 
S(x-x') 



n=l+l 



4ttSk 2 (x)' 



(B8) 

(B9) 
(B10) 

(Bll) 
(B12) 



for K > 0, where k n = n\J~K. In the following, only the case of K < is explicitly presented. The case of K > is 
easily obtained by discretizing the variable k — * k n and by an interpretation of the integrals / °° dk — ► VKJ2^ = i +1 , and 
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by 8(k — k') — ► 5 nn i etc. For our purpose, we need the explicit forms of Xi only up to I — 4: 

X0 = S^)' (B13) 
X\ = —j^=^ - [— Ck{x) sin kx + kSnix) cos kx] , (B14) 

V^Vi(fc) Sjf 0) 

X 2 1 3 {[3- (fc 2 + 2iT)S K 2 (x)] sinfca;- 3fc5 , Ar (a;)CK(x)cosfca;} , (B15) 

\/N 2 (k) S K {x) 

X 3 = 7==^- 4 . - |C x (x) [-15 + 6(fc 2 + A-JSx 2 ^)] sinfcz 
^N 3 (k) (a;) 

+ kS K (x) [15 - (fc 2 + llK)S K 2 (x)] coskx} , (B16) 

I h 05 



X 4 = — { [105 - 15(3fc 2 + 8A")Sk- 2 (jc) + (fc 4 + 35iffc 2 + 24if 2 )S , K 4 (x)] sinfcx 



- kS K (x) [105 - 10(fc 2 + 5K)S K 2 (x)] cos fcz} , (B17) 

where the function 

cosh (xy/^ltj , K < 0, 
<! I- ^ (B18) 



cos (xVk) , K > 0, 



has the properties 



Expanding the density contrast 5(x,9,<f)) in terms of the normal modes, 

U2 



Ck 2 (x)+KS k 2 (x) = 1, (B19) 
dCW(*) = _ x5k(x) . (B20 ) 



<KM,<£) = EE ^s lm (k)x l (k,x)Y l m (e^), (B2i) 

l=0m=-l J ^ 

J; m (fc)=47r y dx5 , K 2 (a;) ^ sin 9d9d<j)5(x, 9, <f>)Xi{k, x)Y l m *{0, 4>), (B22) 
the power spectrum P(k) is defined by 

(s: m (k)S llml (k')) = (27:) 3 S w d mml 6{k ~ 2 k ' ) P(k), (B23) 

where only diagonal elements survive due to the statistical homogeneity and isotropy of the universe. The density contrast 
in redshift space is given by (Matsubara 2000) 

8^{x{z) 1 9 1 ct>)=D{z)\b{z) + f{z) ^ + a(z) A (A + 3K)- 1 } 6(x,6,(t>), (B24) 

where 5(x,9,<p) is the linear density contrast at the present time, b(z) is the linear bias parameter at redshift z, f(z) is 
the logarithmic derivative of the growth factor defined by equation (3.8). The linear operator (A + 3-ftT) -1 denotes the 
Green's function of the operator A + 3K, and 

C K [x(z)} f d\n[D(z)f(z)<l>(z)} \ 
a{z) = SriM] V + dlnS K [x{z)] | ' (B25) 

where ^(z) is the selection function per comoving volume. The selection function $(z) and the redshift distribution 
n(z) = dN(< z)/dz of galaxies in the sample with a fixed sky area are related by n(z)H(z) oc &(z)Sk 2 (x(z)). Therefore, 
the quantity a(z) is more conveniently expressed in terms of n(z): 

a{z) = H(z)± HH(z)D(z)f(z)n(z)]. (B26) 

In this equation, the term Hd\n(HDF)/dz has the order of the inverse of the Hubble scale, therefore is negligible, unless 
the clustering on Hubble scales is calculated. On the other hand, the term Hd\nn/dz is the order of inverse of the scale on 
which the distribution function n(z) varies. When the correlation function on such a scale needs to be obtained, the factor 
a should be retained. Since the selection function <P{z) usually is a decreasing function and the volume factor Sk 2 {x{z)) is 
an increasing function, the distribution n(z) — dN/dz does not significantly vary in the useful redshift range of a sample, 
unless galaxies near redshift edges are included. Therefore, if the redshift distribution n(z) is approximately constant in 
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a sample, the factor a can be neglected. However, when the distribution dN/dz significantly varies in a sample, a must 
be kept in the analysis. 

It is straightforward to obtain the Green's function since the density contrast is expanded by the eigenfunction of 
the Laplacian by equation (B21). The correlation function in rcdshift space is therefore given by ^ s '(zi, Zj, Qij) = 
(5^(x(zi),9i,<pi)5^(x(zj),0j,<pj)), where 9^ is the angle between the two directions (Qi,4>i) and (Qj,4>j)- Defining the 
functions, 



*W( X) = { kHk Xl (k x)P(k) 

* {x) ~ s K 2n -Hx) J 27T* JHfc 2 -4/o n ,( ' ' w ' 



SK ln -\x)J 271-^ k(k' z - AK) n - K ' ' " " ( ' B27 ' ) 
we can explicitly represent the linear formula of the two-point correlation function in rcdshift space (Matsubara 2000): 

2 2n 

^(zuzj^ij) = bibjDiDj^^^ix^xjAj^M, (B28) 

n=0 1=0 

where we abbreviate bi = b(z{), Di = D(zi), Xi = x(zi) etc., D{z) is the linear growth factor normalized by D{z = 0) = 1, 
and x(z) is the comoving distance given by equation (3.4). The quantity X{j is the comoving separation of the two-points, 
which is obtained from x i} Xj, and % by the geometric relation, 

SK 2 (xij) = S K 2 (xi) + S K 2 (xj) - 2C K {xi)CK{x J )SK{xi)S K {x J )cos,e ij 
-KS K 2 {x l )S K 2 {x ] ) (1 + cos 2 %) . 



The coefficients cj"^ are given by 



„(0) 



4 1} = 



1 + \{Pi + ft) + (l + 2 cos 2 %) 

A + Pj + j^PiPj (4 + 3 cos Bi 



M _ 



\K\S K 2 {x i0 ) - -faPjAijAji cos%, 



Aij (cos -fij — 2 cos jji cos Qij ^ + Aji (cos jji — 2 cos jij cos Q^ 

1 



PiAij cos 7^ + PjA^ cos jji 
1 
5 

1 



= & cos 7^ - - + /?,• ( cos Tjj - 



1 



2 4— ~ 
- + - cos 2 % - (cos 2 7^ + cos 2 Tjj) + 4 cos 7^ cos jji cos % 



c (2) - 

.(2) _ 



.(2) _ 



c (2) - 

Co — 



„(2) 



ftft (|jr|S* 2 (a;y) - A^ Aji) \K\S K 2 {x l0 ), 
Pi(3j (A^ cos-fij + Aji cos jji) \K\S K 2 (xij), 



— 29 — 
cos 2 Qij — — + — (cos 2 jij + cos 2 7^) + 10 cos jij cos jji cos 0j 



\K\S K 2 (x v ) 



+ PiPjAijAj, ( cos7^ cos 7^ + - cos% ) , 



Aij (b cos jij cos 2 jji — cos jij + 2 cos jji cos 0y ^ 

+ Aji (j> cos jji cos 2 jij — cos 7j j + 2 cos 7^ cos % ^ 

1 2 — — 
- + - cos 2 Qij — (cos 2 7y + cos 2 jji) + 4 cos jij cos jji cos Qi 

+ 7 cos 2 jij cos 2 jji] , 



where the abbreviation /3j = is employed and /3(z) = f(z)/b(z). The quantity is defined by 



(B29) 

(B30) 
(B31) 

(B32) 

(B33) 

(B34) 
(B35) 

(B36) 

(B37) 

(B38) 
(B39) 



The quantity Aji is similarly defined with the replacement i <-» j. The quantity jij is an angle between the line of sight 
of Xi and the direction of the separation Xij , which can be obtained by the equation 



COS Jij 



S K (xi)C K (xj) - C K (xi)S K (xj) cosOi 



Sx{Xij) 

and jji is similarly defined with the replacement i <-> j. Finally, the quantity Qij is defined by 



cos Qi 



CK(xi)Ci({xj)cosQij + KSK{xi)SK{xj) _ sin 7^ sin 7^ 



C K {xi)C K {xj) + KS K (xi)S K (xj)cos9ij C K (xij) 



- COS Jij COS Jji. 



(B40) 



(B41) 
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The matrix Oij is symmetric, while Aij and 7^ are not symmetric. 

Although the expression of the two-point correlation function (B28) is somewhat tedius, the numerical calculation is 

straightforward. Especially, once the single- variable functions f(z), D(z), x(z), b(z), ^(z), ^ n \x) are calculated and 
tabulated beforehand, then the evaluation of the correlation function does not require any further numerical integrations. 
This property is the essential part for the fast computation of the correlation matrix we have developed in this paper. 

WINDOW FUNCTIONS AND THE EPANECHNIKOV KERNELS 

In this appendix, window functions for a series of kernel functions, i.e., Epanechnikov kernels are explicitly given. 

The Three-dimensional Case 

First we consider the three-dimensional space. The window function is the 3-dimensional Fourier transform of the kernel 
function: 

/p OO 
d 3 xW R (x)e- ikx = 4tt J x 2 dxW R (x)j Q (kx) (CI) 

There are two popular spherical kernels, i.e., the top- hat kernel, 

W PW = ^rs & ( R -^ W^{kR)^-^ 1 (C2) 



where 



is the Heavisidc step function, and the Gaussian kernel, 

wjfHz) = J^J^se-^ 2 ^, W^(kR) = e-^l\ (C4) 

One of the advantages of the top- hat kernel is that it has finite support. It is a disadvantage that the Gaussian kernel 
extends to infinite volume, which should be truncated in actual analyses. However, the Gaussian kernel has an advantage 
that the window function W^(kR) drops off exponentially for large wavenumbers k, and thus the numerical integration 
with this factor is very stable. The function is an oscillating function and the envelope drops off as (kR)~ 2 , which 

converges much slower than the Gaussian window. 

The slow drop-off of the Fourier transform of the top-hat kernel comes from the sharp discontinuity of the kernel at the 
edge, x = R. Accordingly, it is sometimes advantageous to use a kernel that is continuous at the edge, and at the same 
time, has finite support. One of such kernels is the Epanechnikov kernel, 

"Pv-M 1 ™>=w^ « c5 » 

which has a parabolic profile and the envelope of the Fourier transform drops off as (kR)~ 3 . While the Epanechnikov 
kernel is contiuous on the edge, the radial derivative is not. The Epanechnikov kernel can be generalized such that the 
radial derivatives on the edge are also continous. The generalized kernel is called m-weight Epanechnikov kernel, which 
is given by 



(El) 

The original Epanechnikov kernel is the 1-wcight kernel, W R . Top-hat kernel is also from this sequence with m = 0. 
The m-weight Epanechnikov kernel is countinous up to (m — 1) derivatives on the edge. The envelope of the Fourier 
transform drops off as (kR)~ m ~ 2 , so that the convergence is faster for higher orders. The effective width of the kernel is 
smaller than R, especially for higher weights m. In this respect, the size of the kernels is more conveniently represented 
by the "variance" of the kernels, 

1/2 /ON 1/2 



n = 



J d 3 xW R Em \x)x 2 



2m + 5 



R. (C7) 



For the top-hat kernel, the smoothing radius i?T corresponds to Rt = (5/3) 1//2 7?.. In the limit m — > 00, with 1Z fixed, the 
m-weight Epanechnikov kernel reduces to the Gaussian kernel. In fact, 

1Z fixed 

where the Gaussian smoothing length Rq corresponds to Rg = TZ/y/3. Therefore, the generalized Epanechnikov kernels 
contain both the top- hat kernel and the Gaussian kernels at m = and m = 00, respectively. 
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The Two-dimensional Case 
Next wc summarize the case for the two-dimensional space. The window function is given by 

W{W S ) = j ' d 2 6WeMz~ a ' 6 '-Stt J 6deWo s {6)Ml6). (C9) 
The top-hat kernel and its window function are 

= -L^ey. - e), w^(ie s ) = 2 -^, (cio) 

and Gaussian counterparts are 

W M(0) = -^ e - e2 ^ 2 \ W™(W B ) = e~ l ^l\ (Cll) 
2tt6 s 

The Epancchnikov kernel in 2-dimension is given by 

<>(*) = ^(l-J) 0(9 S - 6), = ^« (C12) 

which has a parabolic profile and the envelope of the Fourier transform drops off as (W s )~ 5 / 2 . The m-weight Epanechnikov 
kernel is derived as 

w^\e) = ^(i-^ye { e s -6), ^)(^) = 2m+1(m ( +^r l( ^ ) - (ci3) 

The envelope of the Fourier transform drops off as (fci?)~" l ~ 3//2 , so that the convergence is faster for higher weights. The 
effective width of the kernel is smaller than R, especially for higher weights to. The size of the kernels is alternatively 
represented by the "variance" of the kernels, 

V2 a 



tf s ^ d 2 6W^ m, {6)9 2 = - s (C14) 
J 3 J Vm + 2 

For the top-hat kernel, the smoothing radius 9t corresponds to 6*t = V^^s- In the limit m — > oo, with i) s fixed, the 
to- weight Epanechnikov kernel reduces to the Gaussian kernel. In fact, 

K m) (°) ( C15 ) 

# s fixed 

where the Gaussian smoothing length 8q corresponds to 0q = $s/\/2- 
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-0.1 0.1 -0.15 -0.1 -0.05 0.05 0.1 0.15 

Fig. CI. — The smoothed (R = 15 fo _1 Mpc) correlation function in redshift space in the concordance model: i?M0 = 0.3, J?ko = 0, w = — 1, 
h = 0.7, /baryon = 0.15, b = 1, <Ts = 1. Contour lines indicate the value of the smoothed correlation function around the centers at redshifts 
Zi = 0.04, 0.1, 0.2, 0.3, 0.4, 0.5, 1, 1.5, 2. The coordinate axes correspond to (zj sinSij, Zj cosBij), i.e., the direct 2-space. The contours have the 
intervals of A log 10 ^ = 1/3. Thick contours indicate £jj = 0.1, 0.01, 10 -3 , 10~ 4 from inner to outer contours, respectively. The zero-points 
of £ij are plotted by dotted lines. In the right panels, negative regions are also shown by dashed contours with the same contour levels with 
minus signs. 
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x [h _1 Mpc] x [h _1 Mpc] 

Fig. C2. — The correlation function in real space, t;(x). The thick lines correspond to the concordance model which is the same as in 
the Figure CI. The baryon fraction /baryon is varied in the left panel, and the total density parameter J?mo is varied in the right panel, as 
indicated in the figures. Higher baryon fraction and lower density parameter both give higher amplitudes of the correlations. Dotted lines 
indicate the negative correlations. 
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Fig. C3. — Angular correlation function for the concordance model. The solid line shows the exact prediction of the linear theory (see 
text). The dashed line shows the prediction from Limber's equation. 
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i i i r 
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Fig. C4. — The smoothed angular correlation function u(9) around the smoothing scale (0 B = 10 arcmin.). Thin solid line: the bare angular 
correlations. Thick solid line: smoothed correlations with Limber's equation. Short-dashed line: Taylor approximation up to 1st order. 
Long-dashed line: Taylor approximation up to 2nd order. The lower panel shows the fractional differences with respect to the prediction of 
Limber's equation, which gives the presice value of u{6) in the angle range here. 
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Fig. C5. — A schematic explanation of the projection by the S/N eigenmodes, or KL eigenmodes. Solid ellipses represent the signal 
correlations and dashed ellipses represent the noise correlations in two-dimensional data space. The directions of the S/N eigenmodes are 
indicated by black arrows. The matrix Q\ diagonalizcs the noise correlations, and Q2 rescales the data space to obtain the prewhitened noise 
matrix. Although there is not a unique way to choose a prewhitcning matrix Q, the resulting projection matrix P is independent of the 
choice. The generalized eigenvalue equation directly gives the projection matrix P. The original data space is projected onto the reduced 
data space, retaining only modes with the highest S/N ratio. 



